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Abstract 

In this paper we use the angular deficit scheme [V. Borrelh, F. Cazals, and J.-M. Morvan, Computer Aided 
Geometric Design 20, 319 (2003)] to determine the distribution of Gaussian curvature in developable cones 
(d-cones) [E. Cerda, S. Chaieb, F. Melo, and L. Mahadevan, Nature 401, 46 (1999)] numerically. These 
d-cones are formed by pushing a thin elastic sheet into a circular container. Negative Gaussian curvatures 
are identified at the rim where the sheet touches the container. Around the rim there are two narrow bands 
with positive Gaussian curvatures. The integral of the (negative) Gaussian curvature near the rim is almost 
completely compensated by that of the two adjacent bands. This suggests that the Gauss-Bonnet theorem 
which constrains the integral of Gaussian curvature globally does not explain the spontaneous curvature 
cancellation phenomenon [T. Liang and T. A. Witten, Phys. Rev. E 73, 046604 (2006)]. The locality of the 
compensation seems to increase for decreasing d-cone thickness. The angular deficit scheme also provides 
a new way to confirm the curvature cancellation phenomenon. 

PACS numbers: 46.70.De, 68.55.-a, 46.32.+X 
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I. INTRODUCTION 



Over the past few years, there has been a marked and still increasing interest in the geometrical 
and mechanical properties of developable cones or d-cones, illustrated in Figure [B The d-cone 
offers an especially simple example of focusing of mechanical stress in a thin elastic sheet owing to 
the interplay of bending and stretching energy deformation. Several puzzles regarding the scaling 
of this focusing remain unresolved[l]. Chief among these is the observed vanishing of mean 
curvature at the supporting rim of the d-cone. This paper investigates the possible role of the 
Gauss-Bonnet theorem in explaining this vanishing phenomenon. 




Figure 1: A typical simulated d-cone formed by pushing the center of a hexagonal elastic sheet O against 
a cylindrical container with force F. It has side length / = 60a, container radius R = 38a, deflection 
e = 0.095 and thickness h = 0.102a, where a is the lattice spacing. OP is one boundary of the 60 degree 
sector opposite to the buckled region. Gaussian curvatures are integrated over the shaded region in Sec. III. 
The cylinder is represented by a circular potential shown in Sec. II. 

In principle the shape of a d-cone, including the focusing phenomena noted above^ may be 
found by solving the von Karman equations describing large deformations of thin plates [i2i]. Cerda 
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Figure 2: Top: The radial profile of the d-cone along OP. The slope of the fitted line is -0.092. Bottom: 
The residual plot of difference between the data points and the fitted line. This plot clearly shows that radial 



curvatures are not zero, especially near the rim where y^x^ + = 38. 



and MahadevanyO analyzed the geometry of a single d-cone formed by pushing a thin elastic cir- 
cular sheet into a cylindrical frame by applying a centered transverse force directed along the axis 
of the cylinder, in the limit of pure bending. They obtained the shape of the d-cone by minimizing 
the bending energy and showed that the aperture angle of the buckled part has a universal value 
of 139 degrees and that the sheet exerts a concentrated point force to the frame at the two take-off 
points where the sheet bends away from the frame. This finding was confirmed by Liang and 
Wittenii 

These authors [|5D also uncovered a striking spontaneous curvature cancellation phenomenon as 
noted above. Though the ideal shape of a d-cone has zero radial curvature Crr where r and Q are 
defined as radial and angular components in the material coordinate system, the real shape must 
have nonzero radial curvature at the rim, demonstrated in Figure [2l The localized force requires 
an outward radial curvature of the surface so that the surface elements at the frame may maintain 
equilibrium. As the thickness of the elastic sheet goes to zero, they found that within the numerical 
accuracy, the radial curvature Crr and the azimuthal curvature Cqq are equal and opposite at the 
rim, so that the mean curvature virtually vanishes there. Though it arises from mechanical equilib- 
rium, the phenomenon is purely geometric: it does not involve material parameters [ll|,l5t]. The von 
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Karman equations must be able to account for the cancellation effect, but this has not been done 
yet. 

Any outward radial curvature Crr combined with the inward azimuthal curvature Cqo nec- 
essarily creates a nonzero Gaussian curvature K = Crr Coo. The vanishing of mean curvature 
{Crr + Cqo)/ 2 mcaus that K at the rim approaches a nonzero constant as the thickness of the sheet 
goes to zero. This contrasts with the overall behavior of thin sheets, which approach the isometric, 
developable shape for which it must go to zero. 

The vanishing mean curvature at the rim amounts to a nonlocal geometrical constraint. It 
is nonlocal because it depends on the overall geometry of the sheet. The curvature does not 
vanish when one modifies this geometry by cutting the sheet radially or by replacing the flat 
sheet by a cone[5]. The effect also depends on the nonlocal forcing at the vertex. Given the 
importance of nonlocal geometry, it is natural to look for a connection with the main known 
geometric constraint on curvature: namely the Gauss-Bonnet theorem. This theorem states thatlQ] 
the integral of Gaussian curvature over a region M of a surface is related to the integral of the 
geodesic curvature /^^ over the boundary of that region through j^K dA + Kg ds = 2ti. The 
geodesic curvature for a point on a curve lying on a surface is defined as the curvature of the 
orthogonal projection of the curve onto the tangent plane to the surface at the point. As noted 
below, there is an exact counterpart of the theorem for polyhedral surfaces. 

The Gauss-Bonnet theorem implies that any region of net negative Gaussian curvature, such as 
the d-cone rim, must be compensated by a region of positive curvature elsewhere in the surface or 
a change in the boundary integral. Our aim in this paper is to determine where the rim curvature 
is actually compensated. Since the vanishing mean curvature is the result of global forces in the 
sheet, it would be natural to find that the rim Gaussian curvature is compensated by positive Gaus- 
sian curvature concentrated in the core region. Such nonlocal compensation arises, for example, if 
one pushes a stretched membrane such as a drum head at an interior point. Here positive Gaussian 
curvature appears at the forcing point. Under weak deformation, as one moves away from the 
forcing point, the deformation becomes conical, with no Gaussian curvature. The compensating 
negative Gaussian curvature appears only at the constrained boundary of the membrane-the same 
place where the point force is compensated. It is natural to expect analogous behavior for the 
d-cone. 

These Gauss-Bonnet integrals are expected to be influenced by the elastic thickness h of the 
sheet. As this thickness goes to zero relative to the size R of the container, the sheet approaches the 
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Figure 3: Left: Density plot of the Gaussian curvature K in the d-cone shown in Figure [T] The point P is 
at the lower left corner. The gaussian curvature at each point is indicated by the scale at right. White areas 
indicate off-scale positive or negative values of K. Right: The radial profiles of K as a function of r for 
two different values of thickness h. K is normalized by which is defined in Sec. II. 

isometric state in which the energy penalty for stretching becomes arbitrarily large. In this limit 
the Gaussian curvature must approach zero at almost every point, since Gaussian curvature entails 
stretching [ItI]. Thus the negative contribution to the Gauss-Bonnet integral near the rim must go 
to zero. Although this negative contribution becomes small, it must nevertheless be compensated 
somewhere on the surface. The compensation could occur far from the rim or close to it. We may 
investigate this question by evaluating the integral of K in a band near the rim. If the nagative 
region near the rim is nearly compensated by a positive region within the band, then J^^^^ K <C 
hand 1^1* instead the compensation occurs largely outside the band, then J^^^^^K/ J^^^^ \ K\ 
remains finite as the thickness goes to zero. 

Our numerical study reported below finds that the rim Gaussian curvature is compensated lo- 
cally, not globally. There is a narrow annulus near the rim where the Gauss-Bonnet integral nearly 
vanishes, due to nearly cancelling positive and negative contributions. There is an analogous local 
region around the core. The locality appears to increase for decreasing membrane thickness. Our 
methods allow us to confirm the vanishing-mean-curvature phenomenon in a new way, using the 
Gaussian curvatures. The numerical models we used are presented in Sec. II. In Sec. Ill, we 
describe our numerical findings in detail. In Sec. IV we discuss limitations and implications of 
our findings. 
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II. NUMERICAL METHODS 



A. Numerical Model 

We model an elastic sheet by a triangular lattice of springs of un-stretched length a and spring 
constant k after Seung and Nelsonlsl]. Bending rigidity is introduced by assigning an energy of 
J(l — fii ' 712) to every pair of adjacent triangles with normals hi and 712. When strains are small 
compared to unity and radii of curvature are large compared to the lattice spacing a, this model 
is equivalent to an elastic sheet of thickness h = a^f^JJk made of an isotropic, homogeneous 
material with bending modulus k = J\pij2, Young's modulus Y = 2kajh\pi and Poisson's ratio 
z/ = 1/3. The lattice spacing a is set to be 1. The shape of the sheet in our simulation is a regular 
hexagon of side length R^. The typical value of is 60a, as indicated in Figure [B 

To obtain a single d-cone shape, we need to simulate the constraining container rim and pushing 
force. The rim lies in the x — ^ plane and is described by equation + = i?^, where R is the 
radius of the container. The deflection of the sheet is defined as e = rf/i?, where d is the distance 
by which the center of the sheet has been pushed into the container. A d-cone with deflection e 



touches the container at distance R\J\ + away from the apex. For each d-cone, we define a 
constant as = (e/ (i?vT+~?))^, where ej (i?vT+~?) is the azimuthal curvature Cqq of a 



cone with opening angle 2 tan~^((i/i?) at r = i?vT+~?. In a real d-cone, is slightly different 
from the square of the azimuthal curvature at the rim due to lattice effect and the finite thickness 
h of the sheet. should approach C^q at the rim as h goes to zero. Pushing in the center of the 
sheet is accomplished by introducing a repulsive potential of the form t/force(^i, ^1) = — (^1 + 
a)G{xi,yi)F, where G{xi,yi) = [(1 + + {yi/0'^)]~^^ ^i' Vi zi are the coordinates 

of the lattice point in the center, ^ is a constant of order 0.1a and F is the magnitude of the pushing 
force. G{xi,yi) is introduced to make sure that when the sheet is being pushed the lattice point in 
the center does not stray away from the axis of the cylindrical container, i.e. (xi, yi) (0, 0). The 



constraining rim is implemented by a potential of the form Unm = ^ Cp/[{^/c^^^ — i?)^ + ]^, 
where Cp is a constant and the summation is over all lattice points with coordinates (x^, yi^Zi). The 
value of Cp is choosen so that the shortest distance between the lattice points and the rim is close 
to one lattice spacing. 

The conjugate gradient algorithm [|9[] is used to minimize the total elastic and potential energy 
of the system as a function of the coordinates of all lattice points. This lattice model behaves like 
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a continuum material provided that the curvatures are everywhere sufficiently small compared to 
1/a. This restricts the values of deflection e of d-cone to be below 0.25 [|5|]. e is between 0.09 and 
0.15 in our simulations. 



B. Evaluation of curvatures 

There are two ways to determine the curvatures. 

One way, denoted as Scheme 1, is to obtain the curvature tensor from each triangle in the sheet. 
Once the curvature tensor is known, Gaussian curvature and mean curvature are the determinant 
and average of the diagonal elements of the tensor, respectively. For this scheme to work, we have 
to take the curvature tensor to be constant across each triangle. We calculate it following |iqI] 
using the relative heights of the six vertices of the three triangles that share a side with the given 
triangle. 



Gaussian curvatures can also be obtained by using the angular deficit schemelUy based on the 
Gauss-Bonnet Theorem. The angular deficit scheme approximates the Gaussian curvature of a 
smooth surface S from the angular deficit of an inscribed polyhedral surface. The angular deficit 
of a vertex of a polyhedron is the amount by which the sum of the angles of the faces at the 
vertex falls short of 27r. Let p be a vertex of a polyhedral surface and let its nearest neighbors to 
ht pi,i = 1, . . . , n in clockwise direction. In triangle p PiPi+i, let 7^ denote the angle at p, i.e., 
Zpippi^i , with Pn^i = pi. The angular deficit at p is defined as 27r — Xir=i 7^* equilateral 
triangular lattice model, if a vertex p is not on the edge of the polyhedral surface, then n = 6, and 
the angular deficit and the Gaussian curvature K of the smooth surface S aip satisfy: 



i=i ~ i=i 

in the limit of small lattice spacing [fill]. Since angular deficits can be easily calculated in our 
triangular lattice model, this is a convenient way to evaluate Gaussian curvatures. For a vertex p 
on the edge, n = 3 or 4. The angular deficit at p is ^(tt/S — 7^) where we restrict the angles 
to those subtending the sheet. After this adjustment, the sum of the angular deficits over all the 
vertices on the sheet is equivalent to the sum of angular deficits over all the triangles and should 
be exactly zero. In our simulation, the magnitude of this sum is on the order of 10~^^, which is 
consistent with zero given our numerical round-off errors. 

Equation ([T]) can be seen as the counterpart of the Gauss-Bonnet theorem on smooth surface for 



polyhedral surfaces. In the polyhedral case, the Gaussian curvature of the surface is concentrated 
on the vertices. It is zero on faces and edges. By approximating the polyhedron by a sequence of 
smooth surfaces, one verifies that K consists of 5-function contribution at each vertex j each being 
equal to its angular deficit: K{f) = ^j[^^(7r/3 — 7ij)]5^(r — fj). This special case of Gaussian- 
Bonnet theorem states that the sum of the angular deficits over all the vertices of a polyhedron 
is 27rx, where x is the Euler characteristic of the polyhedron. For example, if a polyhedron is 
homeomorphic to a sphere, then its Euler characteristic x = 2 and the total angular deficits of all 
its vertices is An. So Equation ([T]) simply states that in our regular triangular lattice, the Gaussian 
curvature of a smooth surface can be approximated by the normalized Gaussian curvature of the 
discrete lattice. There is great temptation to generalize this to other kinds of lattices, but it is only 
true when the geometry of the lattices is precisely controlled|ll]. 

These two schemes of obtaining the Gaussian curvature give consistent values everywhere on 
the d-cone except in the core region where the high curvatures invalidate the smoothness assump- 
tion of Scheme 1. All Gaussian curvatures presented in subsequent sections are those of our poly- 
hedral surface determined by the angular deficit scheme, the robustness and accuracy of which has 



been tested extensively 



fly 



III. RESULTS 

Numerically, Figure [3] shows the Gaussian curvature profiles along a radial line OP that is 
57r/6 away from the symmetrical axis of the buckled region as indicated in Figure [B Negative 
Gaussian curvatures are observed around r = 38a in the material coordinate, where the sheet 
touches the confining rim, just as stated in the Introduction. This negative Gaussian curvature 
decays as one moves inward or outward from the rim, but the decay is not monotonic. It first goes 
to zero and then turns positive before approaching zero again. As expected, Gaussian curvatures 
change very dramatically in the core region. 

The existence of two regions with positive Gaussian curvatures around the contact region seems 
to indicate that the negative Gaussian curvatures in the contact region is compensated locally. To 
test this, on the d-cone with thickness h = 0.102a we integrate the Gaussian curvature over a band 
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Figure 4: The integral of Gaussian curvature Jj K within each band Aj as a function of j, where j represents 
the material distance from band Aj to the center of the d-cone. 

within the unbuckled region shown as the shaded region in Figure [B The results are as follows: 

/ / K{e, r)r dr d9 = -5.5 x 10-^ (2) 

Jo=57t/6 Jr=20a 
77r/6 r59a 

/ \K{9,r)\rdrd9 = 2.8 x 10-^ (3) 

_=57r/6 Jr=20a 

These integrals of the Gaussian curvature are two orders of magnitude smaller than the integrals 
of the absolute value of the Gaussian curvature, which means that the negative Gaussian curvature 
at the rim is effectively compensated by the positive Gaussian curvature close to the rim, i.e., the 
Gaussian curvature is compensated locally. 

By comparing the radial profiles of the Gaussian curvature in Figure [3l one would notice that 
the negative Gaussian curvature at the rim decays to zero within a shorter distance on the thinner 
d-cone. Also, on the thinner d-cone the range of two regions of positive Gaussian curvature is 
correspondingly smaller. This suggests that the locality of the Gaussian curvature compensation 
increases as the thickness of the d-cone decreases. 

The local compensation of Gaussian curvature applies not only to the unbuckled region, but 
also to the whole d-cone. We divide a d-cone into circular bands with each band Aj containing 
lattice points that satisfy (i — 1) < r < j, j = 1, . . . , 61, where r is the material distance from 
the lattice point to the center of the d-cone. Then within each band j, we integrate the Gaussian 
curvature and denote the integral as J. K, This ^. K is plotted against j in Figure |4l From Figure|4l 

9 




K/K 




0.05 



0.15 0.2 

h/a 



Figure 5: The normalized average gaussian curvature K at the rim as a function of the thickness of the sheet 
h. Gaussian curvatures at rim points were inferred from interpolation. These interpolated K values were 
then averaged over the contacting rim. 

We can roughly decompose the d-cone to 3 major regions, the core region with < r < 20a, the 
rim region with 20a < r < 50a and the edge region with 50a < r < 60a where the integrals are 
primarily contributed by the buckled region. 
In the core region, we have 



/ / K{9,r)rdrd9 = 2.5 x 10"^ 

J 0=0 Jr=0 

/27T p20a 
/ \K{9,r)\rdrd9 = 1.0 x 10"^ 



(4) 
(5) 



and in the rim region, 



/ / K{9,r)rdrd9 = -1.5 X 10- 

Je=0 Jr=20a 



27r /.50a 



\K{e,r)\rdrde = 2.7 x 10" 



(6) 
(7) 



=0 Jr=20a 

These integrals of the Gaussian curvature are at least one order of magnitude smaller than the 
integrals of the absolute value of the Gaussian curvature in both the core and the rim regions, 
which again suggests that overall the Gaussian curvature is compensated locally. 

Through Gaussian curvatures, we can confirm the vanishing-mean-curvature phenomenon. The 
vanishing of mean curvature {Crr + Coo)/ 2 means that the Gaussian curvature K remains equal 
to Cqq. Since, as stated in Sec. II, Cqq at the rim approaches the constant Kq as the thickness h 
of the sheet goes to zero, the vanishing of the mean curvature would further imply that at the limit 
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h = 0, K = Kq the rim. In Figure [51 there is a clear trend that K/ Kq approaches 1 as /i goes to 
zero. 

IV. DISCUSSION 

If an external normal point force is applied to an elastic membrane, positive Gaussian curvature 
is necessarily created at the forcing point. In general this must be compensated by a region of 
negative Gaussian curvature elsewhere in the membrane, owing to the Gauss-Bonnet theorem. 
For the case of the d-cone lattice considered here the content of the theorem is quite precise: the 
angular deficits defined above must sum to zero. Thus any deficit due to an external force must 
be exactly compensated elsewhere in the sheet. This compensation may occur close to the applied 
force or far from it, depending on the case. The size of the compensating region can depend on 
the magnitude of the applied force, the forces at distant boundaries and the elastic thickness of the 
membrane. Similarly the degree of compensation falls off with distance from the source in a way 
that depends on the case. The nature of compensation for weak point forces is tractable using the 
linearized von Karman equations. 

Our interest here is the case of the d-cone, where puzzling nonlinear and nonlocal effects like 
the vanishing of mean curvature at the supporting rim are seen. As explained above, this vanishing 
appears to require nonlocal interaction between the rim region and the rest of the system. Any 
nonlocal compensation of Gaussian curvature would provide a potential mechanism for this non- 
local interaction. However, our numerical analysis shows that the rim gaussian curvature is almost 
entirely compensated close to the rim. Indeed, the compensation appears to become more local as 
the elastic thickness of the membrane is reduced. Similar locality appears in the compensation for 
the central point force. Thus we have found nothing to suggest a coupling of the central region to 
the rim via the Gauss-Bonnet theorem. 

Our findings here are modest and limited. We have only explored one structure: the d-cone, 
although the locality of compensation is of general interest and relevance. Clearly many variants 
of the d-cone could be studied, as well as many other cases of point forces on curved and strained 
membranes. We expect that the increased locality we observed with decreasing thickness obeys 
some form of scaling behavior, but we have not attempted to determine it numerically or analyt- 
ically. In order to demonstrate such scaling, any numerical measurement would need to have a 
finite element scheme with much more dynamic range than the one used here. 
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The puzzle of vanishing mean curvature at the rim of a d-cone is as open as ever after this study. 
Our work in preparation investigates the generahty of this curvature cancellation and explores 
simplifying assumptions to enable understanding. 
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